Application of modified simplex method to biomagnetic inverse problem 
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1 Introduction 

Difficulty in determination of source brain activities 
for observed biomagnetic signals (magneto¬ 
encephalograms, MEGs) is known as the 
biomagnetic inverse problem. Moving Dipole 
Method [1] has been developed for solution of this 
problem. Two difficulties have been pointed out for 
this approach. The first problem is that the number of 
sources must be known before the analysis. However, 
it is difficult to fix the number of sources for the 
brain activities which change from time to time. The 
second problem arises for the case of a relatively 
large source. The estimated location of the source 
tended to deviate from that of the real source [2]. We 
have developed a new algorithm, called Moving 
Mesh Method (MMM) [3][4] in order to overcome 
these problems. 

One of the difficulties in MMM was that it required a 
large number of iterative calculations. It was also 
difficult to find good initial conditions for iterative 
calculations. The first problem was solved by 
introducing the modified Simplex Method (SM) [5] 
where the dipoles were assumed at vertexes of 
tetrahedrons which can be transformed so that it 
gives the best goodness of fit (GOF) (see Flowchart). 
We solved the second problem by introducing 
Singularity Search Method (SSM) [6] to locate 
dipoles responsible for the radial components of 
magnetic fields observed on imaginary spheres as 
singularity points of an indicator function given by a 
product of the sphere radius and the radial magnetic 
component. At the same time, the number n of 
dipoles was assumed by comparing the moment of 
estimated dipoles (see phantom experiments). 

The feasibility of the new algorithm was tested in 
phantom experiments where magnetic signals were 
generated by current electrodes and the estimates of 
magnetic sources determined by MMM for the 
magnetic signals were compared with the real ones, 
and also in MEG experiments where somatosensory 
evoked fields (SEFs), the auditory evoked fields 
(AEFs) and the epileptic fields were analyzed. 


2 Methods 

2.1 Flowchart 

The brain was assumed as a homogenous sphere, and 
the magnetic fields generated by dipoles were given 
as [7]: 

B=L(r) P . (1) 

Where B : a set of the magnetic field 
L : the lead field matrix 
r : the dipole positions 
P : the dipole moments. 

For given dipoles, we can estimate P using the 
inverse matrix of L. L may include errors if the 
assumed dipoles deviate from the real dipoles. The 
errors can be minimized, by re-estimating the P after 
repositioning the assumed dipole and iterating this 
procedure, until the difference between the observed 
and predicted magnetic fields from the estimated 
dipoles becomes minimal. A flow chart of the 
method is shown below, and is described in details in 
next section. 

Stepl: Estimate r and n to account for the radial 
component of B observed on imaginary 
spherical surfaces by SSM. 

Step2: Estimate P using the regularized g-inverse 
matrix (equation 2) for the regularized 
parameter (equation 3) determined by the 
Generalized Cross Validation (GCV) 
method [8]. 

Step3: Determine GOF given by the dipoles 
estimated in step 2. 

Step4: Move a dipole according to SM algorithm 
and recalculate P and GOF. 

Step5: Finalize the solution for P and r for those 
giving the maximum GOF. 

2.2 Estimation of dipole moment 

In Step 1 one needs to estimate r and n for B, which 
is a typical inverse problem contained in the 
electromagnetism. This problem can be simplified if 
the solution is constrained to a problem to know r 
and n for the radial component of B observed on 
imaginary spherical surfaces. The solutions are 
found as singularity points of an indicator function 
given by a product of the sphere radius and the radial 
magnetic component. 



Step 2 to determine P is ill-posed due to the ill-posed 
nature of L. We solved this problem by applying the 
technique of the Tikhonov regularization [9]. And 
the moments are determined as: 

P=(L t L+A G t G) _1 L t B . (2) 

Where A : regularization parameter 
G : penalty matrix. 

G is a penalty matrix to avoid solutions closer to the 
sensors [10]. A was determined by GCV method, so 
that V(A) becomes minimal, which is given by 

V(X)=|{I-(L+L’)C}B| 2 

/ {trace{I-(L+L’)C}} 2 (3) 

where C=(L t L+XG) _1 L t 

L’: the difference matrix between the lead 
field of a dipole and that of the nearest one. 
Simulation study indicates that V(A) is a unimodal 
function, and therefore, the Golden Section method, 
known as the fast algorithm, was used to determine 
the A in equation 3. 

2.3 Estimation of dipole position 

Since B is a nonlinear function of r, r has to be 
estimated for B by an iterative algorithm. Step 1 
estimating dipoles for the radial components of B 
gave initial conditions for this iterative estimation of 
r. The estimate was iteratively improved using the 
modified SM algorithm, assuming dipoles on 
vertexes of tetrahedrons which were reflected, 
expanded or contracted to maximize P. 

2.4 Recording system 

Biomagnetic recording was conducted using the 
Shimadzu Biomagnetic Imaging System constructed 
of 201 sensors arranged over the whole head which 
is shown in Figure 1. This system introduced the 
Flux Locked Loop • FLL • circuits using the flux 
modulation by the alternating current (AC) bias. 



Figure 1: The Shimadzu Biomagnetic Imaging 
System. 


3 Results 

3.1 Phantom experiments 

The precision of the new algorithm was studied in 
phantom experiments with single and two dipoles 
producing the magnetic fields with intensities 
comparable with the normal biomagnetic fields. Data 
were sampled at lKHz, band-passed (1-100 Hz) and 
averaged for 100 trials. The initial and final results of 
the iterative estimates of the dipoles are shown in 
Tables 1 and 2 in comparison with the real dipoles. 
In the single dipole case SSM revealed a single 
major tetrahedrons and four minor tetrahedrons 
whose moments were much smaller (<1/10) than the 
major one. Improvement of the estimates by SM 
algorithm was conducted on the major tetrahedron. 
Similarly, in two dipoles case, improvement by SM 
algorithm was conducted on two major tetrahedrons 
revealed by SSM. Altogether, four and eight dipoles 
at vertexes of one and two tetrahedrons were 
estimated for a single and two dipole cases, 
respectively. They were distanced so close each 
other (the maximum distance, 0.1 mm). The errors of 
the dipole estimates (the distance between the 
gravity of the estimated dipoles and the real dipoles) 
were less than about 5 mm, indicating that the new 
method is capable of correctly estimating the 

Table 1: The estimated result from the data 
obtained from the phantom with single dipole. SL 
represents the positional error of the estimated from 
the real dipole. _ 



X 

Y 

Z 

Px 

Py 

Pz 

5L 

Initial 1 

4.5 

49.0 

35.1 

5.5 

15.8 

-17.8 

4.8 

Initial 2 

-43.4 

83.9 

18.4 

-0.5 

-0.1 

-0.7 


Initial 3 

-43.3 

84.3 

18.8 

0.5 

0.1 

0.7 


Initial 4 

18.5 

96.4 

-19.0 

0.5 

0.1 

1.0 


Initial 5 

17.0 

99.3 

-20.1 

-0.9 

0.0 

-0.6 


Final 1 

3.3 

52.9 

38.4 

0.9 

2.8 

-3.7 

3.2 

Real 1 

1.2 

52.3 

36.1 

3.8 

14.2 

-19.6 



Table 2: The estimated result from the data 
obtained from the phantom with double dipoles. 



X 

Y 

Z 

Px 

Py 

Pz 

§L 

Initial 1 

-0.5 

-27.8 

-52.0 

-20.9 

-18.0 

12.3 

6.2 

Initial 2 

48.0 

-15.5 

-35.3 

-10.9 

-8.7 

-10.1 

17.0 

Initial 3 

49.2 

-40.4 

-58.3 

1.1 

2.5 

-1.1 


Initial 4 

-40.8 

-74.1 

-16.7 

1.8 

-0.9 

-0.3 


Initial 5 

-18.5 

-85.1 

-22.7 

-1.7 

0.2 

0.7 


Final 1 

2.6 

-22.8 

-58.4 

-5.4 

-5.2 

1.9 

5.1 

Final 2 

29.9 

-14.1 

-30.5 

-2.9 

-2.5 

-3.6 

2.0 

True 1 

3.2 

-27.7 

-57.0 

-7.4 

-14.1 

6.3 


True 2 

32.1 

-14.7 

-29.4 

-32.9 

-13.7 

-23.3 



































multiple dipoles with no priori knowledge of the 
number of dipoles. 

3.2 Application to biomagnetic fields 

The feasibility of the new method was further tested 
in three types of MEG experiments. Informed 
consent was obtained from all subjects prior to the 
study. The first case was the sensory-evoked fields 
(SEFs) recorded from a faculty of the laboratory 
(male, aged 39). Electrical stimuli were delivered to 
the left median nerve located at the wrist (duration: 
0.3 msec, strength: 10 mA), MEGs were recorded 
using a bandwidth filter (1-100 Hz), and averaged 
for 100 stimulation trials (the acquisition window 
was limited to 256 msec with stimulation at t=50 
msec). 

MMM was conducted on the data recorded at t=72 
msec when the brain activity was maximal. The 
results are shown in Figure 2. In this case, four 
dipoles on a single tetrahedron were estimated, 
clustering in the contralateral (right) primary 
somatosensory cortex to the stimulated hand. 

The second case was the auditory evoked fields 
(AEFs) recorded from the same subject. Tone bursts 
(frequency, 1000 Hz) were delivered to the subject’s 
left ear. MEG recording was conducted in the same 
way as for the SEFs. MMM conducted t=120 msec 
(Figure 3) revealed twenty dipoles which were 
clustered in two major groups, i.e. one at the 
contralateral (right) primary auditory cortex to the 
stimulated ear, and the other at the ipsilateral (left) 
cortex. 

The last case was the epileptic fields from a patient 
(female, aged 52). Magnetic resonance imaging 
(MRI) of the patient’s brain revealed a lesion in the 
right hippocampus. MEGs and electroencephalogram 
(EEG) were recorded simultaneously (sampling rate, 
500 Hz and band-pass filter, 1-100 Hz). MMM 
conducted at the time of the epileptic spike four 
dipoles clustered around the hippocampal lesion. 

4 Discussion 

Phantom experiments showed that MMM was 
capable of correctly estimating current sources of the 
model brain signals with no prior knowledge of the 
current sources. MMM also revealed major current 
sources for SEFs and AEFs in a restricted focus of 
the contralateral primary somatosensory and 
audiotory cortical areas, respectively. MMM also 
was capable of localizing the current sources for the 
epileptic fields which were noisier than SEF and 
AEF. Epileptic sources were scattered in a much 
larger area than for SEF and AEF. This may 
represent scatter distribution of epileptic sources. 



Figure 2: Estimated dipoles from the response to the 
stimulation of the left median nerve. The vector 
arrows’ direction and length represent the dipole 
direction and magnitude respectively. The slice oj 
the MRI was selected as that at which the estimated 
source had the largest moment. 


a) 



Figure 3: Estimated dipoles from the response to the 
stimulation of the left ear. a) The slice of MRI was 
selected as that at which the estimated source in the 
right hemisphere had the largest moment, b) The 
slice of MRI was selected as that at which the 
estimated source in the left hemisphere had the 
largest moment. 



Figure 4: Estimated dipoles from the epileptic spike. 
The slice of the MRI was selected at the right 
hippocampus. 










Importantly, the computation time for the new 
algorithm of MMM was 1/10 (5 min in case of SEF 
in Figure 2 using a DEC workstation model 600au). 

Acknowledgements 

We thank Dr. T. Kobayashi and Dr. Y. Kida, 
Department of neurosurgery, Komaki City Hospital, 
for recording the epileptic data. 

This work was supported by Special Coordination 
Funds for Promoting Science and Technology of 
Japan. 

References 

1. B.N. Cuffin, “A comparison of moving dipole 
inverse solutions using EEG’s and MEG’s”, 
IEEE Transactions on Biomedical Engineering 
11, 905-910, 1985. 

2. R. Hari, “On brain's magnetic responses to 
sensory stimuli”, Journal of Clinical 
Neurophysiology , 8, 157-169, 1991. 

3. S. Kajihara, T. Tomita, S. Okamura, 

A. Arakawa, T. Tomita, Y. Yoshida, 

H. Ohshima and Y. Takanashi, “Spread source 
analysis by moving mesh method”, in Recent 
Advances in Biomagnetism , T. Yoshimoto, M. 
Kotani, S. Kuriki, H. Karibe, and N. Nakasato, 


Eds. Sendai: Tohoku University Press, 1999, pp. 
236-239. 

4. S. Kajihara, T. Tomita, Y. Kondo, A. Arakawa, 

S. Okamura, T. Tomita, Y. Yoshida and 
Y. Takanashi, “Moving mesh method for 
reconstructing some spread sources in the brain”, 
Brain Topography , in press. 

5. J.A. Nelder and R. Mead, “ A simplex method 
for function minimization”, Computer journal 7, 
308, 1965. 

6. K. Yamatani, and K. Ohnaka, “ An estimation 
method for point sources of multidimensional 
diffusionequation”, Appl. Math. Modelling , 21- 
2, 77-84(1997) 

7. J. Sarvas, “Basic mathematical and electro¬ 
magnetic concepts of the biomagnetic inverse 
problem ”, Phys.Med.Biol., 32:11-22, 1987. 

8. G.H. Golub, M. Heath and G. Wahba, 
“Generalized cross-validation as a method for 
choosing a good ridge parameter”, 
Technometrics 21-2, 215-223, 1979. 

9. A. N. Tikhonov, “Solution of incorrectly 
formulated problems and regularization method”, 
Soviet Math. 4, 1035-1038, 1963. 

10. K. Sekihara, B. Scholz, and H. Bruder, 
“Separating of correlated and uncorrelated 
source activities”, Advances in Biomagnetism 
Abstracts , 233-234, 1993. 



